      subroutine nmrvf(ln,ecmn,ucentr,cntfug,gridx,nx,
     >   jdouble,njdouble,s1,s2,reg,jstart,jstop)
      dimension jdouble(njdouble),ucentr(nx),reg(nx),gridx(nx),
     >   cntfug(nx)

      x = gridx(jstart)
      dx= gridx(jstart+1) - x
      f2 = ucentr(jstart)+cntfug(jstart)-ecmn
      h1=dx*dx
      h2=h1/12.0

      if (jstart.eq.1) then
C  We get here if the integration is going to start from the first X point.
C  This means that S1 is the solution of the differential equation at X=0
C  S2 is the solution at the first X point. 
         s1=0.0
         t1=0.0
C  LN = 1 is a special case
         if (ln.eq.1) t1=-h1/18.0
         if (ecmn.ne.0.0) t1=t1*ecmn
      else
         j=jstart-1
         f1=ucentr(j)+cntfug(j)-ecmn
         t1=s1*(1.0-h2*f1)
      end if

      reg(jstart) = s2
      t2=s2*(1.0-h2*f2)
      
      istart=2
      do while (jstart.gt.jdouble(istart).and.istart.lt.njdouble)
         istart=istart+1
      end do
      istart=istart-1
      istop=njdouble-1
C  JDOUBLE(ISTART) points to the first doubling of DX that happens after JSTART
C  JDOUBLE(ISTOP) points to the last doubling of DX that happens before JSTOP
      
C    integration loop
      do i=istart,istop
         j1=max(jstart,jdouble(i))+1
         j2=min(jstop,jdouble(i+1))
         do j=j1,j2
            f3 = ucentr(j)+cntfug(j)-ecmn
            t3 = 2.*t2-t1+h1*f2*s2
            s3 = t3/(1.0-h2*f3)
            reg(j)=s3      
      
            t1=t2
            t2=t3
            f0=f1
            f1=f2
            f2=f3
            s0=s1
            s1=s2
            s2=s3
         end do
         dx=2.0*dx
         h1=dx*dx
         h2=h1/12.0
         t2=s3*(1.0-h2*f3)
         t1=s0*(1.0-h2*f0)
      end do
      return
      end

